AD-A164  1*5 
UNCLASSIFIED 


UNSTEflDV  SOLUTION  OF  THE  BOUNDRRV  LAVER  EQUATIONS  HITH  1/2 
APPLICATION  TO  SPA.  .  (U)  AIR  FORCE  INST  OF  TECH 
HRIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENGI. .  K  J  LANGE 
JUN  83  AFIT/GAE/ENV/84D-11  F/G  28/13  NL 


AFIT/GAE/ENY/84D-1 1 


0 


DTIC 


UNSTEADY  SOLUTION  OF  THE  BOUNDARY 
LAYER  EQUATIONS  WITH  APPLICATION  TO 
SPACE  SHUTTLE  TILES 
THESIS 

Karen  J.  Lange 
Captain,  USAF 

AFIT/GAE/ENY/9  4D-1 1 


Approved  for  public  release;  distribution  unlimited 


AFIT/GAE/ENY/84D-1 1 


UNSTEADY  SOLUTION  OF  THE  BOUNDARY 
LAYER  EQUATIONS  WITH  APPLICATION  TO 
SPACE  SHUTTLE  TILES 

THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Aeronautical  Engineering 


Karen  J.  Lange,  B.S. 
Captain,  USAF 


June  1985 


Accesion  For 


NTIS  CRA&I 
OTIC  TAB 
Unannounced 
Justification 


□ 

□ 


By . 

Di^t  itution/ 


Avaitab&ty  Codes 


EJist 


m 


Avail  and/or 
Sp&da# 


Approved  for  public  release;  distribution  unlimited 


ACKNOWLEDGMENTS 


I  would  like  to  recognize  the  contributions  of  the  following 
individuals,  whose  guidance  and  support  helped  me  complete  this  project. 

First,  I  want  to  thank  Maj.  James  K.  Hodge  for  his  long  hours  of 
assistance  and  understanding  during  the  many  weeks  that  I  was  struggling 
to  get  positive  results  from  my  computer  program. 

Secondly,  I  wish  to  express  my  appreciation  to  Dr.  James  Hitchcock, 
who  motivated  me  to  pursue  heat  transfer  by  the  quality  of  his  guidance 
and  his  wisdom  in  the  classroom.  His  classes  were  very  interesting  and 
I  enjoyed  them  immensely. 

Last,  but  not  least,  I  want  to  thank  my  husband,  Joel  and  my  dog, 
Beau,  for  their  many  long  months  of  devotion  and  understanding  without 
which  I  could  not  have  finished  this  project. 


< 


( 


iii 


< 


V,  «*.  • 


>V%.V-V>V'  V-V-  O 

i  1  M-&,  M-h  fcAV-V  b 


Contents 

Acknowledgments  . 

List  of  Figures  . 

List  of  Tables  . 

List  of  Symbols  . 

Abstract  . 

I.  Introduction  . 

Background  . 

Objectives  and  Approach  . 

II.  Theory  . 

Governing  Equations  . 

Low  Temperature  Viscosity  . 

Heat  Transfer  Coefficient  . 

Oblique  Shock  Theory  ........ 

III.  Finite  Difference  Solution  . 

Grid . 

Numerical  Method  . 

Momentum  and  Energy  Equations  .  .  . 

Metrics  . 

Optimized  Successive  Over-Relaxation 

Continuity  Equation  . 

Heat  Rate  . 

Boundary  Conditions  . 

Wall  Boundary  Conditions  . 

Edge  Boundary  Conditions  . 

Initial  Conditions  . 

IV.  Results  . 

Program  Verification  . 

Non  Isothermal  Cases  . 

Transient  Cases  . 

Pitching  Plate  Case  . 

V.  Conclusions  and  Recommendations  .  .  . 

Conclusions  . 

Recommendations  . 

Appendix  A:  Dorodnitsyn  Transformation  .  .  . 

Appendix  B:  Thesis  Computer  Code  . 

Bibliography  . 


List  of  Figures 


igure 

1  Thin  Skin,  HRSI,  or  FRSI  Panels  in  Wedge  Test  Article  .  .  . 

2  Illustration  of  Y-Value  Determination  for  the  Optimized 

Blasius  Grid  . 

3  n  versus  Velocity  in  the  Transformed  Plane  . 

4a  Fixed  Grid . . . 

4b  Expanded  View  of  Fixed  Grid  . . 

5a  Physical  Grid  at  Steady  State  . 

5b  Expanded  View  of  the  Physical  Grid  at  Steady  State  . 

6  Comparison  of  Blasius  Solution  with  the  Numerical  Solution  . 

7  Thin  Skin  Wedge  Heating  Data  Comparison  with  Predictions 

for  Isothermal  Wall  . 

8  Nonisothermal  Wall  Solution  at  14  Degrees  with  a  HRSI  Wall 

Temperature  of  800  Degrees  . 

9  Nonisothermal  Wall  Solution  at  10  Degrees  with  a  HRSI  Wall 

Temperature  of  800  Degrees  . 

0  Nonisothermal  Wall  Solution  at  3  Degrees  with  a  HRSI  Wall 
Tanperature  of  800  Degrees  . 

1  Nonisothermal  Wall  Solution  at  14  Degrees  with  a  FRSI  Wall 

Temperature  of  700  Degrees  . 

2  Transient  Nonisothermal  Wall  Solutions  at  14  Degrees  .... 

3  Comparison  of  Transient  Nonisothermal  Solutions  with  the  Wall 

Temperature  a  Function  of  the  Heating  Rate  . 

4  Deflection  Angle  Time  History  . 

5  Comparison  of  Nonisothermal  Pitching  Plate  Solutions  for 

14  Degrees,  9  Degrees,  and  3  Degrees  . 

6  Comparison  of  the  Wall  Temperatures  for  the  Nonisothermal 
Pitching  Plate  at  14  Degrees,  9  Degrees,  and  3  Degrees  .  .  . 

7  Vhre£  Time  History  for  X/L  Location  of  .016667  . 


60 


.  ”«  ■ 


18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 


Wall  Temperature  Time  History  for  the  Streamwise  Location 


of  .016667  . 

h/href  Versus  Deflection  Angle  for  X/L  Location  of  .01667 

h/href  Time  History  for  X/L  Location  of  .466667  . 

Wall  Temperature  Time  History  for  X/L  Location  of  .466667 
h/href  Versus  Deflection  Angle  for  X/L  Location  of  .466667 

h/href  Time  History  for  X/L  Location  of  .483333  . 

Wall  Temperature  Time  History  for  X/L  Location  of  .483333 
h/href  Versus  Deflection  Angle  for  X/L  Location  of  .483333 

hA*ref  Time  History  for  X/L  Location  of  .566667  . 

Wall  Temperature  Time  History  for  X/L  Location  of  .566667  . 
h/href  Versus  Deflection  Angle  for  X/L  Location  of  .566667 

h/href  Time  History  for  X/L  Location  of  .983333  . 

Wall  Temperature  Time  History  for  X/L  Location  of  .983333  . 
h/href  Versus  Deflection  Angle  for  X/L  Location  of  .983333 


.61 

.63 


.64 


.65 


66 


.67 


.68 


69 


.70 


.71 


73 


.74 


.75 


vi 


,-V7 


List  of  Symbols 


a  Speed  of  sound 

Cp  Specific  heat  capacity  at  constant  pressure 

H  Total  enthalpy 
h  Heat  transfer  coefficient 

hiso  Heat  transfer  coefficient  for  an  isothermal  wall 
hre£  Reference  heat  transfer  coefficient 
hx  Local  heat  transfer  coefficient 
h  g  Heat  transfer  coefficient  with  .9  recovery  factor 
I  Number  of  spatial  node  points  in  X  direction 
i  Node  point  in  X  direction 

J  Number  of  spatial  node  points  in  y  direction 

j  Node  point  in  y  direction 

k  Conductivity 

L  Plate  length 

M  Mach  number 

M  Freestream  Mach  number 
00 

Nux  Local  Nusselt  number 
P  Pressure 

PQ  Stagnation  pressure 
Pr  Prandtl  number 

q  Heat  rate 

qQ  Heat  rate  at  the  wall 

R  Gas  constant 

r  Recovery  factor 


RC  Time  constant 


Re  Reynolds  number 

Rex  Local  Reynolds  number 

T  Temperature 

Taw  Adiabatic  Wall  temperature 

Te  Boundary  layer  edge  temperature 

Tw  Wall  temperature 

TQ  Stagnation  temperature 

T.  Freestream  temperature 

T*  Reference  temperature 

t  Time 

U  Dorodnitsyn' s  transformation  of  streamwise  flow  velocity 
u  Streamwise  flow  velocity 

V  Dorodnitsyn' s  transformation  of  normal  flow  velocity 

v  Normal  flow  velocity 

X  Dorodnitsyn' s  transformation  of  streamwise  spatial  coordinate 
x  Streamwise  spatial  coordinate 

Y  Dorodnitsyn' s  transformation  of  normal  spatial  coordinate 

y  Normal  spatial  coordinate 

a  Deflection  angle  of  experimental  plate 
3  Wave  angle 

6  Boundary  layer  thickness 

y  Specific  heat  ratio 

e  Wedge  angle 

n  Transformed  Y  in  the  computational  plane 
u  Viscosity 

ue  Boundary  layer  edge  viscosity 


ix 


Freest ream  Viscosity 
p  Density 

Pe  Boundary  layer  edge  density 

po  Freest ream  density 

£  Transformed  X  in  the  computational  plane 
e  Streamwise  location  of  wall  temperature  step 
Superscripts 
n  Time  level 

s  Iteration  level 

*  Reference  temperature  conditions 

'  Dimensional  Quantity 

Subscripts 

e  Boundary  layer  edge  condition 


i,j  Spatial  node  point 
°°  Freestream  condition 


Abstract 


Wind  tunnel  tests  were  conducted  on  Space  Shuttle  Orbiter 
insulating  articles.  The  data  was  processed  by  a  heat  estimation 
program,,  revealing  a  large  discrepancy  in  magnitude  of  the  heat  transfer 
coefficient  over  the  Space  Shuttle  tile  when  compared  to  that  yielded  by 
flat  plate  theory  and  thin  skin  test  results.  An  initial  suggestion  for 
the  cause  of  this  discrepancy  was  the  nonisothermal  wall  effect. 
Understanding  this  effect  was  the  main  objective  of  this  project. 

In  order  to  investigate  these  effects,  an  unsteady  compressible 
boundary  layer  code,  THESIS,  was  developed.  First  the  boundary  layer 
equations  were  transformed  using  the  Dorodnitsyn  transformation  to  allow 
the  equations  to  be  solved  similar  to  incompressible  equations.  Then, 
the  equations  were  then  transformed  to  the  computational  plane  and 
approximated  by  finite  differences.  Optimized  successive  over¬ 
relaxation  is  then  used  to  solve  the  difference  equations.  The  grid 
used  for  this  program  is  optimized  based  on  the  Blasius  solution,  but 
the  program  is  not  limited  to  this  because  it  has  a  generalized  grid. 

Upon  running  this  program,  leading  edge  startup  problems  arose. 
Therefore,  a  similarity  solution  is  assumed  at  the  leading  edge  which 
corrected  the  problem. 

The  program  was  verified  with  an  incompressible  case  which  was  in 
good  agreement  with  the  Blasius  solution,  and  a  compressible  case  with 
an  isothermal  wall  which  agreed  well  with  approximate  theory.  The 
isothermal  case  showed  the  solution  was  dependent  on  the  grid.  The 
steady  state  and  transient  nonisothermal  wall  cases  displayed  the  same 
trend  as  the  approximate  nonisothermal  solution.  All  cases  had  the  same 


dramatic  decrease  in  the  heat  transfer  coefficient  and  the  subsequent 
recovery  in  the  streamwise  direction.  The  most  significant  results  came 
from  pitching  the  plate  with  a  nonisothermal  wall.  A  large  lag  due  to 
the  nonisothermal  wall  effects  exists.  These  results  indicate  that  the 
nonisothermal  wall  effects  may  be  affecting  the  wind  tunnel  results. 


UNSTEADY  SOLUTION  OF  THE  BOUNDARY 
LAYER  EQUATIONS  WITH  APPLICATION  TO 
SPACE  SHUTTLE  TILES 


I.  INTRODUCTION 


The  Space  Shuttle  has  completed  numerous  missions  in  a  relatively 
short  time  and  has  become  an  important  tool  in  the  study  and 
understanding  of  space.  The  unique  feature  of  the  shuttle  is  its 
reuseability.  The  Thermal  Protection  System  (TPS)  protects  the  shuttle 
from  extreme  temperatures  upon  reentry  into  the  atmosphere  and  enables 
it  to  be  reused.  Understanding  how  the  TPS  responds  to  the  temperature 
flucuations  is  important  for  future  applications  such  as  designing  new 
hypersonic  vehicles.  Studies  can  be  accomplished  through  many  means 
such  as  actual  flights,  wind  tunnels,  and  computers. 


Background 

Flight  test  maneuvers  were  simulated  in  a  twenty  inch  hypersonic 
wind  tunnel  to  verify  a  technique  used  to  estimate  heating  rate  and 
thermal  parameters  from  Space  Shuttle  thermocouple  data.  The  tests  were 
conducted  by  the  Air  Force  Wright  Aeronautical  Laboratory  at  Wright- 
Patterson  AFB,  Ohio.  Three  types  of  materials,  a  thin  skin  stainless 


steel  plate,  and  two  TPS  panels,  a  high-temperature  reuseable  surface 
insulation  (HRSI),  and  a  flexible  reuseable  surface  insulation  (FRSI), 
were  flush  mounted  on  a  test  plate  with  a  wedge-shape  leading  edge  as 
shown  in  Figure  1.  The  wedge  was  1.25  feet  by  7  inches  wide,  and  a  6 


1 


inch  by  six  inch  panel  interfaced  with  the  stainless  steel  surface  at  a 
streamwise  distance  of  7  inches  from  the  leading  edge  or  at  x/L  =  .467. 
Different  thermocouple  instrumentation  was  installed  in  each  panel.  The 
deflection  angle  (o$  was  held  at  a  constant  angle  or  varied  temporally 
between  3  degrees  and  14  degrees  during  a  wind  tunnel  run.  The  wedge 
was  pitched  at  a  constant  pitch  rate  of  2  degrees/second  to  selected 
deflection  angles  and  held  for  a  prescribed  time. 

The  conditions  in  the  test  section  were  known  from  the  test  section 
Mach  number  of  14.237,  the  stagnation  temperature  of  2020  degrees  R,  and 
the  stagnation  pressure  of  1600  pounds  per  square  inch.  The  air  was 
assumed  to  be  a  perfect  gas  with  a  specific  heat  ratio  (y  )  of  1.4. 

Fifty  nine  tests  were  run.  The  first  26  were  with  a  HRSI  test 
article,  27-45  were  with  a  FRSI  test  article,  and  the  last  14  test  runs 
were  with  the  thin  skin  article  which  simulated  an  isothermal  wall.  The 
thin  skin  results  indicate  a  good  agreement  with  the  isothermal  wall 
predictions  but  the  HRSI  and  FRSI  panels  differed  markedly  from  Eckert 
flat  plate  theory.  The  main  difference  between  the  TPS  panels  and  the 
thin  skin  stainless  steel  panel  was  the  thermal  properties.  The  TPS 
panels  had  a  much  lower  conductivity  so  the  surface  temperature 
increased  much  more  rapidly  than  the  stainless  steel.  Therefore,  the 
tile  was  much  hotter  than  the  stainless  steel  wedge,  causing  a 
temperature  discontinuity.  This  nonisothermal  wall  was  indicated  by  the 
results  as  the  cause  of  an  additional  thermocouple  lag  during  the 
pitching  maneuvers.  (7:5-7)  It  also  causes  a  large  decrease  or  increase 
in  laminar  heating  depending  on  the  direction  of  the  temperature  jump. 


The  results  of  the  wind  tunnel  tests  show  that  nonisothermal 


effects  should  not  be  neglected  for  vehicles  with  laminar  flow  that  have 
numerous  interfaces  of  different  low  heat  capacity  materials. 
Understanding  the  effects  will  be  important  to  understanding  the 
thermocouple  data  from  the  Space  Shuttle  and  in  future  hypersonic 
vehicle  design. 

Objectives  and  Approach 

Hie  nonisothermal  wall  is  caused  by  the  transient  response  of  the 
wall  temperature  which  is  caused  by  the  limited  duration  test  and  by  the 
motion  of  the  plate.  Understanding  these  nonisothermal  wall  effects  and 
the  heating  rate  is  important  to  understanding  how  the  TPS  responds. 

One  method  used  to  accomplish  this  is  wind  tunnel  tests.  Another  method 
is  a  two-dimensional  computer  program  using  the  transient  compressible 
boundary  layer  equations.  Ideally,  the  computer  program  would  include 
the  solution  of  the  flow  through  a  shock  wave  and  the  determination  of 
the  heat  rate  with  both  convection  and  conduction  included.  Ultimately, 
it  would  be  able  to  match  a  full  wind  tunnel  profile  so  the  numerical 
results  could  be  compared  to  the  wind  tunnel  results. 

A  two-dimensional  computer  program  with  the  transient  heat  transfer 
boundary  layer  equations  is  written  using  an  implicit  finite  difference 
technique.  A  pitching  flat  plate  will  model  the  wind  tunnel  wedge. 

Using  the  same  free  stream  conditions  as  the  wind  tunnel,  shock  theory 
is  used  to  determine  the  flow  field  behind  the  shock  wave.  The  first 
checkout  of  the  computer  program  is  an  incompressible  steady  state  flat 
plate  with  an  isothermal  wall.  The  results  should  match  the  Blasius 
profile.  To  further  verify  the  program,  an  isothermal  wall  is  run  to 


steady  state  using  the  same  conditions  as  the  thin  skin  stainless  steel 
panel  in  the  wind  tunnel  tests.  This  is  completed  at  deflection  angles 
of  3  degrees,  10  degrees,  and  14  degrees.  The  heat  transfer  coefficient 
ratio  (h/href)  should  compare  well  with  Eckert's  flat  plate  theory.  To 
investigate  the  nonisothermal  wall  case,  temperature  increases  of  170 
degrees  R  for  the  FRSI  panel  and  270  degrees  R  for  the  HRSI  panel  are 
placed  at  x/L  =  .467  with  a  deflection  angle  (a)  of  14  degrees.  The 
results  should  indicate  a  large  initial  decrease  in  the  heat  transfer 
coefficient  ratio.  Moving  down  the  plate,  the  ratio  should  increase  and 
level  off.  To  further  study  the  nonisothermal  effects,  transient  cases 
are  run. 

There  are  two  types  of  transient  cases.  The  first  case  is  at  a 
constant  deflection  angle  of  14  degrees  with  the  wall  temperature 
changing  with  time  by 

Tw"  =  (1  -  e~t/RC)ATw  +  .^-1  (1) 

where  Tw°  and  Tw0-1  are  the  wall  temperature  at  the  present  and  previous 
time,  respectively,  ATw  is  the  temperature  step,  RC  is  the  time 
constant,  and  t  is  the  time.  The  results  should  converge  on  the  same 
solution  as  the  steady  state  nonisothermal  wall  case  explained  above. 

The  second  transient  case  pitches  the  plate  from  14  degrees  to  3  degrees 
and  back  to  14  degrees  with  both  isothermal  and  nonisothermal  walls. 
These  cases  will  show  the  capability  of  the  program  and  will  enhance  the 
knowledge  of  nonisothermal  wall  effects. 


II.  Theory 

Governing  Equations 

The  unsteady  compressible  boundary  layer  equations  were  solved 
numerically  for  both  an  isothermal  and  a  nonisothermal  wall. 

The  Navier-Stokes  equations  for  laminar  flow  are  simplified  by 
assuming  that  gradients  in  the  streamwise  direction  are  small  in 
comparison  with  normal  gradients.  The  normal  velocity  component  is  also 
small.  Terms  of  order  d/L  compared  to  unity  have  been  ignored.  The 
resulting  boundary  layer  equations  (8:17-38)  in  conservative  form  are: 

continuity  3p*  +  3p*  u*  +  3  p 1  v*  =  0  (2) 

3t'  3x'  3y' 

x-momentum  3p' u*  +  3p' u2'  +  Sp'v'u1  =  -  3P*  +  3  / p1 3u' \  (3) 

3t'  3x'  3yf  3x'  3y'  l  3y'  ) 

y-momentum  3£^_  =0  (4) 

3y' 

energy  3p  'H'  +  3p  'u'H  +  3p  *  v'H1  =  3P1  -  3q'  +  3  (  u'  p1  3u'  j  (5) 

3t'  sx'  3y'  3t'  3y'  3y'  V  3y' / 


and  are  subject  to  boundary  conditions  at  the  edge  of  the  boundary  layer 
and  at  the  wall.  The  equations  are  nondimensionalized  using  these 
factors 
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In  order  to  solve  the  unsteady  compressible  boundary  layer 
equations  in  a  form  similar  to  the  incompressible  boundary  layer 
equations,  the  Dorodnitsyn  transformation  (9:2-4)  is  used.  This 
transformation  is 

y 

X  =  x  Y=/Pdy  t=t  H  =  H 

o 

U  =  u  V=pv+^Y  +  U3Y 

at  ax 


where  Y  is  scaled  with  the  density.  The  nondimensional  governing 
equations  become 
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The  equations  look  like  the  incompressible  equations  with  the  density 


imbedded  in  the  viscous  terms  and  the  pressure  term.  For  detailed 
information  on  transforming  the  boundary  layer  equations,  see  Appendix 


The  density  { p')  is  calculated  from  the  perfect  gas  equation  of 


state 


P'  =  p'R'T' 


where  R  is  the  gas  constant  and  P'  and  T'  are  the  pressure  and 
temperature  respectively.  Since  the  boundary  layer  is  assumed  to  be 
thin,  the  pressure  at  the  edge  of  the  boundary  layer  is  transmitted 
without  change  to  the  surface.  Therefore,  the  pressure  is  known 
everywhere  in  the  boundary  layer.  Also,  the  temperature  is  known  from 
the  solution  of  the  energy  equation.  The  only  unknown  is  density  which 
is  then  solved  by 


p'  =  P' 
RT1 


which  nondimens ionally  is 


p  =  r  £ 
y-1  T 


Low  Temperature  Viscosity 

The  standard  relationship  used  to  calculate  viscosity  (  u)  is  that 
given  by  Sutherland's  Law 


u'  =  2.27  x  10“8  T8^'  lb  sec  (21) 

T'  +  198.6  ft2 

According  to  Fiore  (5),  this  may  not  be  valid  for  hypersonic  wind  tunnel 
applications.  In  hypersonic  tunnels,  the  freestream  temperature 
generally  falls  in  a  range  between  =  30  degrees  R  and  =  200 
degrees  R  which  is  outside  the  demonstrated  range  of  applicability  of 
Sutherland's  law.  An  alternate  method  to  calculate  viscosity  is 

u*  =  2.32  x  10-8  T1/2'  lb  sec  (22) 

1  +  220  ft2 

T'x109/t' 

This  equation  provides  good  correlation  to  experimental  data  at  low 
temperatures  and  corresponds  well  with  viscosity  calculations  by 
Sutherland's  Law  at  temperatures  above  T.  =  300  degrees  R  (2:13-18). 

Heat  Transfer  Coefficient 

Heat  transfer  over  a  flat  plate  can  be  calculated  by  Eckert  theory, 
assuming  that  the  flow  medium  is  an  ideal,  single  component  gas  with  no 
dissociation  and  has  constant  properties,  independent  of  temperature  and 
pressure.  Also,  the  surface  temperature  is  assumed  constant.  Using 
these  assumptions,  the  Nusselt  Number  for  laminar  flow  over  a  flat  plate 
(8:134-139)  is 


Nux  =  0.332  Pr1/3  Rex1/2 


(23) 


Pr  is  the  Prandtl  Number  and  is  defined  as 


where  cp  is  the  specific  heat  at  constant  pressure  of  the  fluid,  k  is 
the  conductivity,  and  y  is  the  fluid  viscosity.  Rex  is  the  local 
Reynold's  Number  defined  by 


(24) 


where  the  subscript  e  denotes  boundary  layer  edge  conditions.  Using 
this  definition  of  Nusselt  number, 


(25) 


the  heat  transfer  coefficient  for  an  isothermal  wall  on  a  flat  wedge 
(hiso)  for  laminar  flow  is  given  by 


This  h|SQ  is  independent  of  wall  temperature,  but  the  properties  are  not 
constant.  The  viscosity  and  density  vary  with  temperature  across  a 
boundary  layer.  Therefore,  Eckert's  reference  temperature  (4:10-13) 


T*  =  T'e  +  *5 (T'w  -  T'e)  +  .22(T'aw  -  T’e) 
is  used  as  an  approximation  and 


(27) 


(28) 


where  the  *  indicates  the  condition  at  the  reference  temperature. 


The  heat  transfer  coefficient  can  be  calculated  numerically  at  the 


wall  by 

h.9  *  qi</<-9T,0-  -  T'w) 


where  qQ  is  the  heat  rate  and  is  equal  to 


^o 


-k*  3T'  \ 

3y'  /  y=0 


(29) 


(30) 


The  temperature  gradient  at  the  wall  is  determined  from  the  solution  to 
the  boundary  layer  equations.  The  adiabatic  wall  recovery  factor  is 
always  assumed  to  be  .9  even  though  it  usually  is  defined  as  /Pr.  Using 
/Pr  as  the  recovery  factor  will  only  change  the  magnitude  of  the  heat 
transfer  coefficient  by  the  factor  .9T'0  -  T'w  which  will  usually  be 

^T'o  “  T'w 

small  for  high  speed  flow. 

A  ratio  of  the  heat  transfer  coefficient,  h  g  to  a  reference  heat 
transfer  coefficient  at  zero  deflection  angle,  hre£,  is  determined  where 
the  reference  heat  transfer  coefficient  is 


(31) 


The  conditions  in  front  of  the  shock  wave  are  considered  the  freestream 
conditions  and  are  constant  at  any  deflection  angle.  Thus,  hg  is 
nondimensionalized  by  a  constant  which  is  scaled  by  /x. 

For  a  nonisothermal  wall  on  a  flat  wedge  with  a  wall  temperature 
jimp  from  T^'  to  T2'  at  t,  the  heat  transfer  coefficient  for  laminar 
flow  is  given  by  (8:134-139) 


h.9  =  hiso 


(T'aw-T'i)  +  (l-_«3/4r1/3  (T 


’l-T'2) 


/(.9T'  -T'  )  (32) 


Oblique  Shock  Theory 

The  edge  conditions  of  the  boundary  layer  change  as  the  plate  is 
pitched  from  3  degrees  to  14  degrees.  It  is  necessary  to  know  these 
edge  conditions;  hence,  they  are  calculated  using  oblique  shock  wave 
theory  assuming  adiabatic  flow  of  a  perfect  gas  (11:84-88).  The  first 
task  is  to  calculate  the  Mach  number  after  the  shock,  M2.  In  order  to 
do  this,  the  shock  wave  angle,  8,  is  calculated  by  iteration  using 


tan  e  =  2  cot  8 


M2  sin2  8  -  1 
1 


M2  (y+cos  28)  +  2 


(33) 


where  Mj^  is  the  freestream  Mach  number,  y  is  the  specific  heat  ratio, 
and  a  is  the  wedge  angle.  Once  8  is  known,  the  Mach  number  is 
calculated  using 


f 

r  jczi  2  2  1 

1+2  sin^  8 

-1 

m2  = 

YM2  sin2  8  -  2 

L. 

1/2 


> 


sin2  (8-9) 


(34) 


Using  this  Mach  number,  the  temperature  and  pressure  behind  the  shock 
wave  are  determined  by  the  following  equations 


T'  =  1  +  y— 1  M, 


(35) 


where  T'  is  the  stagnation  temperature. 

The  velocity  is  also  known  because  the  Mach  number  and  the  speed 
sound,  a',  is  known  as 

a'  =/y  RT' 1 

where  R  is  the  gas  constant  for  air. 


III.  Finite  Difference  Solution 


A  two  dimensional  computer  program  is  developed  using  the  transient 
compressible  boundary  layer  equations.  The  equations  are  transformed 
first  with  a  Dorodnitsyn  transformation  and  then  transformed  into  a 
computational  plane.  The  convective  terms  are  differenced  using  second 
order  upwind  differencing  and  the  viscous  terms  are  second  order 
central  differenced.  The  time  terms  use  a  first  order  backward 
difference.  The  differencing  is  done  implicitly  and  successive  over¬ 
relaxation  is  used  to  solve  the  momentum  and  energy  equations.  The 
continuity  equation  is  solved  by  integration  using  the  trapezoidal  rule. 
In  the  program,  a  flat  plate  simulates  the  wind  tunnel  wedge. 

Grid 

The  grid  generation  is  an  important  aspect  of  setting  up  any  finite 

difference  solution.  In  this  case,  the  computational  domain  is  1.25 

feet  in  the  x'-direction  and  approximately  10  times  the  boundary  layer 

thickness  (  6  )  in  the  y'-direction.  Both  x'  and  y'  are 

nondimensionalized  by  the  length  of  the  plate,  L,  so  x  goes  from  0  to  1 

and  y  goes  from  0  to  approximately  106  .  The  computational  domain  has 

L 

1830  grid  points  or  61  points  in  the  streamwise  direction  and  30  points 
in  the  normal  direction.  These  numbers  can  be  varied  to  specify  a  finer 
or  coarser  grid. 

An  optimized  grid  based  on  a  Blasius  solution  at  a  Reynold's  number 
of  2X10^  is  used  as  the  basis  of  this  program's  grid.  It  was  optimized 
to  minimize  truncation  error  in  velocity  gradients.  To  optimize  a  grid 
based  on  a  Blasius  solution,  the  velocity  profile  is  determined  and 
plotted  for  a  specific  location.  Since  similarity  of  velocity  profiles 
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exist  in  an  incompressible  boundary  layer,  any  location  can  be  chosen, 

in  this  case,  x  equal  to  1.0.  A  constant  au  of  .1  was  chosen.  At  each 

increment  of  .1,  the  y  value  was  determined  as  illustrated  in  Figure  2. 

Outside  the  boundary  layer  another  method  which  minimizes  the  change  in 

step  size  is  used  to  choose  y  because  no  Au  exists.  The  y  values  used 

in  the  program  are  given  in  Table  I.  When  transformed  to  the 

computational  grid,  the  plot  of  n  versus  u  would  be  a  diagonal  line  as 

shown  in  Figure  3.  Therefore,  the  truncation  terms,  u  ,  u _ ,  u  ,  H  „ 

3  nn  nnrr  n  nn 

and  Hnnri  would  all  equal  0  in  the  boundary  layer  for  Pr=l. 

This  optimized  Blasius  grid  is  modified  for  a  Reynold's  number  of 
9x10^  which  matches  the  wind  tunnel  condition.  It  is  then  scaled  by  the 
✓x  and  divided  by  the  Prandtl  number  to  allow  for  the  thickness  of  the 
thermal  layer.  In  the  streamwise  direction,  the  node  points  are  equally 
spaced.  Every  X  location  is  specified  in  the  domain  and,  therefore,  the 
X  value  does  not  change  in  the  Y  direction.  These  X,  Y  pair 
combinations  specify  an  initial  grid  in  the  Dorodnitsyn  plane  which  is 
fixed.  Since  this  grid  is  based  on  an  incompressible  solution,  it  must 
be  scaled  to  account  for  the  change  in  the  boundary  layer  due  to 
compressibility  by  a  factor.  A  factor  of  1.0  was  used  for  the  grid 
transformation  in  this  project.  This  is  now  the  fixed  grid.  Since 
density  changes  with  both  temperature  and  time  the  grid  in  the  physical 
plane  also  changes  due  to  the  Dorodnitsyn  transformation  equation  from 
the  fixed  grid  to  the  physical  grid. 

Y 

y  =  /  1  dY  (38) 

o  p 
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Table  I 


Comparing  Fig  4  and  Fig  4b,  the  fixed  grid,  and  Fig  5  and  Fig  5b,  the 
physical  grid  at  steady  state,  the  density  effect  can  be  seen. 


Using  the  fixed  grid,  a  transformation  can  be  made  to  a  rectangular 
computational  grid.  This  grid  is  in  the  £-n  plane  where  n  is  in  the 
streamwise  direction  and  n  is  perpendicular  to  ?.  The  metrics 
n^,  and  Ty  appear  in  the  differential  equations  and  can  be  solved  from 
the  expressions 

5X  =  VJ 

=  -Xn/J 

nX  =  "VJ  (39) 

ny  =  X5/J 

J  =  Vn  - 


Since  X  is  not  changing  in  the  n  direction,  X^  is  0.  This  eliminates 
all  Cy  terms  in  the  differential  equations  and  the  other  metrics 
become 


?x  =  1A5 

Hy  =  lAn  (40) 


The  terms  X^  and  are  determined  using  finite  differences.  The  is 
determined  analytically  which  is  discussed  in  the  section  on  metrics. 

A  unique  feature  of  this  program  is  that  the  grid  used  can  be 
easily  replaced  by  another  grid.  The  derivatives  are  calculated  in  the 
program  so  the  only  replacement  is  the  generation  of  the  fixed  grid. 
This  is  based  on  the  assumption  that  X  does  not  change  as  you  move  in 
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Figure  4a  Fixed  Grid 


19 


EXPANDED  VIEW 


Figure  4b  Expanded  View  of  Fixed  Grid 
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the  Y  direction.  If  X  does  change  in  the  Y  direction,  then  Xn  is  not 
zero,  and  additional  terms  would  have  to  be  added  to  the  differential 
equations.  If  Xn  is  not  zero,  a  variation  from  the  normal  boundary 
layer  assunptions  would  have  to  be  implemented. 

Numerical  Method 


Hie  nondimens ional  transient  two-dimensional  boundary  layer 
equations  are  solved  nianerically  by  using  an  implicit  finite  difference 
scheme  because  it  is  unconditionally  stable  according  to  linear  theory. 
First  the  Dorodnitsyn  transformation  transforms  the  governing  equations, 
Eqs  12  through  15,  and  then  they  are  transformed  from  the  X  -  Y  plane  to 
the  k  -  n  computational  plane.  The  equations  become 


Continuity  cxu^  +  nxUn  +  HyVn  =  0 


(41) 


X-Momentum 


Ut  +  U(U5CX  +  Unnx)  +  VUnny  =  1 


Re 


Mp(UnnY> 


ny 


(42) 


Energy 


Ht  +  U(H55X+Hnnx)  +  VH^  =  1  Pfc  + 

P 


jip  Hnny 


RePr 


HnV|  'V  + 

n 


up(1 


-Pr)  n*  (u/| 


2Re 


(43) 


The  momentum  and  energy  equationsare  solved  for  U  and  H  respectively 
using  SOR  and  then  the  continuity  equation  is  solved  to  give  V. 

Momentum  and  Energy  Equations.  The  momentum  and  energy  equations 
must  first  be  finite  differenced  and  some  assumptions  must  be  made.  On 
a  stationary  plate,  the  pressure  gradients  are  zero  since  inviscid 
oblique  shock  theory  is  assumed.  In  the  momentum  equation,  Eq  42,  the 
dP/dX  term  has  been  assumed  to  be  zero  for  all  cases.  For  a  pitching 


plate,  the  transients  in  the  hypersonic  flow  or  supersonic  edge  flow 
have  a  much  shorter  time  constant  than  the  thermal  response  of  the  TPS 
surface  and  can  be  neglected  (7:5).  For  example,  the  HRSI  tile  has  a 
time  constant  of  three  seconds  and  the  transient  time  constant  is 
.00026  sec.  Also,  no  curvature  is  being  assumed.  The  equations  would 
need  to  be  rederived  to  account  for  curvature. 

Once  the  equations  are  transformed  to  the  5  -  n  plane,  the  first 
order  derivatives  in  the  momentum  and  energy  equations,  Eqs  42  and  43, 
are  differenced  using  three  point  upwind  differencing.  This 
differencing  scheme  guarantees  diagonal  dominance  and  is  second  order 
accurate  in  the  computational  plane.  The  time  derivative  is  differenced 
using  two  point  backward  which  implies  the  scheme  is  only  first  order 
accurate  in  time.  The  second  derivatives,  such  as  the  viscous  terms 


Unn,  Hnn  ,  and  ,  are  approximated  with  second  order  accurate  central 
differencing.  The  like  terms  are  then  combined.  The  equations  are 


n  n 

written  for  U^j  and  H^j  and  become 


Momentum 


n  /  n- 

Uij  =  Uij 


+  25„UAtUi_1j  -  .5  Cj^UAtUi_2 j  + 


2nxUAt  +  2nyVAt  +  nyAt  ( ( pu)  j-n/2  '’y) 


Re 


ij-1 


n  I 


(44) 


[.5nxUAt  +  ,5nyVAt]  U^_2  +  n^At  (  (pu)  j+1/2  V  uij+l  j  /UC0EF 


Re 


UCOEF  =  1  +  35x  UAt  +  3nxUAt  +  3nyVAt  +  rVAt( Duj+i/2ny+Puj-i/2  ny^ 


Energy 


+ 


n  (  n-1  n 

Hij  =|Hij  +  2^X  UAt  Hi-lj  '  *55x  Hi-2j 


[2nxUAt  +  2nyVAt  +  n^t  ( (pp)  j.^V  ] 


RePr 


n  n 

[.SHjjUAt  +  .SHyVAt]  Hi  j_2  +  HyAt  ((p»)  j+1/2V  + 

RePr 


l  dP  +  (l-Pr)ny  At 
p 


2  Re 


2n  2n 

(pu>  j+l/2nYuij+l  ~  ((Pu)  j+l/2rV+(Pll)  j— l/2nY)Ui j  + 


2n 


(py)j_l/2  nyui j-1 


/UCOEF  (45) 


UCOEF  =  1  +  3£x  UAt  +  3nxUAt  +  3nyVAt  +  nyAt(pUj+1/2ny+PMj_1/2  n^)  (46) 


RePr 


Where  i  and  j  define  the  location  in  the  £-n  plane  and  n  is  the  time 
level.  The  notation  for  PPj+1/2  and  puj-l/2  defines  an  average  value  of 
density  time  viscosity  where 

puj+l/2  =  <Pwj+l  +  puj)/2  (47) 

pwj-l/2  =  (ppj  +  ppj-l)/2 

The  momentum  equation  is  a  nonlinear  equation  and  the  energy 
equation  depends  on  the  velocity  calculated  by  the  momentum  equation. 

To  linearize  the  equations  for  the  first  iteration,  the  coefficients 
could  be  lagged  "in  time"  or  "in  space."  The  coefficients  in  this 
program  are  based  on  the  solution  for  U,  V,  and  H  at  the  previous  time 
and  same  location  for  the  first  iteration,  and  subsequently  on  values 


from  the  previous  iteration.  Also,  the  temperature  value  in  the  energy 
equation  is  lagged  in  time  for  the  first  iteration  due  to  its  dependence 
on  H  and  U. 

Metrics.  Imbedded  in  the  coefficients  are  the  metrics  5^,  nx,  and 

n^.  These  quantities  can  be  calculated  using  Eq  40.  A  first  order 

central  difference  is  used  to  determine  X^.  Since  the  points  in  the 

streamwise  direction  are  equally  spaced,  X^  is  a  constant  which  is  equal 

is  and  is  therefore  also  a 

constant.  Two  differencing  schemes  are  used  to  find  Yy  the  inverse  of 

Oy,  depending  on  the  Uy  location  in  the  equation.  For  Hy  inside  the 

viscous  terms,  such  as  |>pny  u  ]  and  [u£ny  H^]  ,  the  Y^  is  differenced 

Re  n  Re  n 

the  same  direction  as  the  pu  term.  For  instance,  the  term  is  a  first 
order  forward  difference  for  Pyj+i/2  anc*  a  f*rst  order  backward 
difference  for  When  the  complete  term  is  differenced,  the 

overall  accuracy  is  second  order.  In  all  other  places,  Y  is  solved 
using  a  second  order  central  difference. 

The  derivative  Y^  needs  to  be  calculated  to  determine  n^.  It  is 
calculated  analytically.  Since  the  grid  is  scaled  by  the  Blasius 
transformation  in  the  X-direction, 


to  the  delta  x.  The  inverse  of  X^ 


nb  =  YiW'  (48) 

M  * 

is  used,  along  with  the  fact  that  x  is  equal  to  £  times  dx,  to  calculate 
Y^.  This  implies  that 


Y  =  rWSdX' 


(49) 


(50) 


dY  =  1  nb/ax  H 

2  - 

/He  /T 

Now,  replacing  nb  with  Eq  48  and  canceling  like  terms 

Yr  =  j_5Y_  (51) 

5 

where  5=i-l.  Therefore,  V  becomes 

Y  r  =  sx  (52) 

(i-1) 

If  a  different  grid  is  used,  another  method  for  calculating  would 
have  to  be  used.  The  method  used  should  be  developed  with  care. 

Initially  Y^  was  calculated  using  a  central  difference  with  a 
backward  difference  at  the  end  of  the  domain.  A  large  amount  of  leading 
edge  error  was  induced  with  this  method.  For  instance,  at  i=2,  the 
analytical  ^  is  .000099437,  and  the  central  difference  is  .000200479 
for  a  dX  of  .016667.  Even  for  a  smaller  dX,  error  is  induced. 

Downstream  of  the  leading  edge,  the  central  difference  is  accurate. 
Therefore,  the  grid  could  be  evaluated  numerically  downstream  and 
analytically  near  the  leading  edge  of  the  flat  plate. 

Once  the  momentum  and  energy  equations  are  linearized  and  the 


metrics  calculated,  the  equations  are  written  as 
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These  equations  are  programmed.  At  the  i=2  and/or  j=2,  a  two  point 
backward  difference  is  used  on  the  convective  terms  if  necessary.  Then 
an  optimized  successive  over-relaxation  technique  is  used. 

Optimized  Successive  Over-Relaxation.  Optimized  successive  over¬ 
relaxation  (SOR)  is  a  technique  which  is  used  in  an  attempt  to 
accelerate  the  iterative  procedure  (1:132-133).  The  arbitrary 
corrections  to  the  intermediate  values  of  the  unknowns  are 


s+1’ 

s+1 

s' 

uij 

=  cu  Uij  + 

(l-o))  Uij 

(57) 

s+1' 

s+1 

s' 

Hij 

-  *  Hij  + 

(l-<0)  Hij 

s+1  s’ 

where  U^  is  the  most  recently  calculated  value  of  U^j,  U^  is  the 

value  from  the  previous  iteration  adjusted  by  the  application  of 

s+1' 

equation  57.  The  term  Ujj  is  the  newly  adjusted  or  "better  guess"  for 

s+1 ' 

Ujj.  The  formula  is  applied  immediately  at  each  point  after  U^j  is 
s+1  s+1 

obtained  and  U-j  replaces  U^.  The  same  iterpretation  is  used  for  the 
H  equation.  In  the  computer  program,  U  is  iterated  once  and  then  H  is 
iterated  once.  This  is  repeated  until  U  and  H  are  converged. 

Omega  ( uj  is  the  relaxation  parameter  and  is  between  zero  and  two. 
Outside  the  boundary  layer,  <o  is  one.  An  optimized  u,  w0pt,  is 
calculated  for  the  locations  inside  the  boundary  layer.  The  u>0pt  is  not 
easily  defineable  for  a  three  point  backward  scheme  but  the  two  point 
backward  scheme  is  defined.  This  is  used  in  helping  to  define  u>0pt 
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the  three  point  backward  scheme  (6).  Since  the  boundary  conditions  are 
not  taken  into  account,  <u0pt  is  not  a  true  optimum.  The  U  velocity 
equation  is  put  in  the  form  of 

AUf j_i  +  BU| j  +  CU^j+^  =  U^j_2  +  ui_ij  +  U^_2j  +  All  other  terms  (58) 

Then  a  maximum  eigenvalue,  x^,  for  Jacobi  iteration  is  determined 
locally  using  the  A,  B,  and  C  coefficients  where 

X|,  =  -2  /NT  cos  ir  (59) 

max  B  JMAX  +1 


and,  therefore 


Ak  =  “2  //2U ijAtnx  +  2Atvij  +  At  Pwj-l/2  \  /  At  Puj+l/2  \  cos  *  <60) 


max 


Y,  ReYn  (Yj-Yj.i)  J  lReYn  (Yj+1-Yj)J 


JMAX+1 


1+l.SUjjAt  +  1.5Ujj  n^At  +  l,5VjjAt  + At  /P^j+l/2  + 


YnRe 


Y^xi-Y 


j+l_xj  Yj“Yj~l 


As  an  approximation,  the  two  circled  2's  are  replaced  with  l's  (6).  The 
At  can  be  divided  out  and  the  x^  become 


xk - 2//2uijnX  +  vij  +  1  pyj-l/2 

max  '  J  J 


1  ph 


j+1/2 


COS  IT 


(61) 


Yn  j?  Re  (Yj-Yj^)!  \YnRe  (Yj+rYj)y 


JMAX+1 


1  +  1.5Ujj  +  1 .5UJ j  nx  +  l.SVij  +  1 


At 


rpu 


j+1/2  +  ppj-l/2 


YnRe  \ Y 


j+l’Yj  Yi-Yj-1 


and 


opt 


1  +  /  l-x2 
k 


(62) 


31 


For  a  real  Jacobi  eigenvalue,  the  A.^  is  limited  to  be  approximately  less 
than  .995,  i.e.  uopt  is  approximately  less  than  1.85.  The  <i)Qpt  is 
recalculated  in  the  program  after  a  specified  number  of  time  iterations. 
The  wQpt  should  be  used  in  the  program  because  it  helps  reduce 
computation  time.  In  some  problems,  it  can  reduce  the  computation  time 
by  a  factor  of  30  which  is  significant  (1:132-133). 

Once  the  updated  values  of  U  and  H  are  calculated  for  every  j 
location,  the  V  velocity  can  be  calculated  from  the  continuity  equation. 

Continuity  Equation.  The  continuity  equation  is  used  to  calculate 
the  V  velocity.  First,  the  equation  is  transformed  to  the  £-n  plane, 

Eq  41,  and  the  V  is  solved  for 


V  =  -  Y  U,  +  Yr  li 
n  n  5  5  n 


Now,  the  equation  is  integrated  to  find  V  using  the  trapezoidal  rule 


n  n 


V-  •  -  v- ■ 
13  13- 


ij-i  =  ~Yn  (ul  ij  +  ij-iy  Y^n  ij  +  un  ij-1^ 
Xj-  \  2  /  Xr\  2  j 


The  are  finite  differenced  with  a  second-order  three  point  upwind 


difference.  The  Un  +  Un  looks  like  Un 

ij  i j— 1  i j— 1/2 

2 


and  is  differenced 


with  a  central  difference  about  j-1/2.  Like  terms  are  combined,  and  the 

differenced  equation  becomes 

n  n  n  n  n 

vij  =  vij-l  +  /~3  Yn  +  YC\ uij  +  Yn  Ui-lj  '  Yn  ui-2j  + 


V  4  X5  X5/ 
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(65) 


/-3  Yn  -  Vc'vUij.!  +  Yn  Ui.ij.!  -  Y„ 


Once  V  is  calculated,  the  iteration  is  complete.  The  U,  H,  and  V 
equations  are  iterated  until  the  convergence  criteria  is  met.  Upon 
convergence,  the  heat  rate  at  the  wall  can  be  calculated  for  that  i- 
location. 

Heat  Rate.  The  heat  rate  at  the  wall,  qQ,  is  defined  by  Eq  30. 

This  equation  must  be  tranformed  to  the  Dorodnitsyn  plane  before  it  can 
be  differenced.  After  it  is  transformed,  the  equation  becomes 


^-X’o'wall 


_3T'  \ 

3Y'  J  Y=c 


where  k  is  the  conductivity  and  is  calculated  by 


k'  =  ii'Cp' 
Pr 


The  derivative  3T/3Y  is  transformed  to  the  e-n  plane  and  a  three  point 
windward  difference  is  used  for  second  order  accuracy  in  the  S-n  plane. 
Since  q  is  a  dimensioned  quantity,  the  T  and  Y  must  be  dimensioned  and 
the  equation  becomes 


l 

•*'  P'wall  Y“3Til  +  4Ti2  “  Ti3^ 
Pr  L'  AY  \  2  / 


where 


AY  =  (-3YU  +  4Yi2  -  Yi3)  /  2 
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Now  the  heat  transfer  coefficient  is  calculated  by 

h.9  '  -  Tiall>  <70' 

where  .9  is  the  recovery  factor  instead  of  /Pr.  Now  this  heat  transfer 

coefficient  can  be  divided  by  the  reference  heat  transfer  coefficient  as 

defined  by  Eq  31. 

Upon  completion  of  this  step,  the  program  can  march  in  the  i 
direction  and  the  solution  procedure  is  repeated  at  this  new  location. 
Finally,  when  the  end  of  the  plate  is  reached,  the  program  goes  to  the 
beginning  of  the  flat  plate  and  starts  again. 

Boundary  Conditions 

In  order  to  completely  define  this  problem,  the  initial  condition 
and  the  wall  and  edge  boundary  conditions  must  be  specified. 

Wall-Boundary  Conditions.  The  wall  boundary  conditions  are  no  slip 
conditions  which  means  U(i,l)  =  o  and  V(i,l)  =  o.  The  wall  temperature 
depends  on  the  case  being  run.  For  an  isothermal  case,  the  wall 
temperature  is  constant  at  530  degrees  R.  Several  nonisothermal  cases 
with  a  temperature  step  at  x/L  of  .466667  were  run.  Again,  this 
temperature  step  depends  on  the  case.  Some  of  the  cases  had  a 
temperature  step  of  170  degrees  R,  increasing  the  wall  temperature  from 
530  degrees  R  to  700  degrees  R  to  simulate  the  FRSI  wind  tunnel  tests, 
and  others  had  a  temperature  step  of  270  degrees  R,  increasing  the  wall 
temperature  to  800  degrees  R  to  simulate  the  HRSI  wind  tunnel  tests.  In 
addition  to  these  nonisothermal  cases,  the  wall  temperature  was  varied 
in  time  from  530  degrees  R  to  800  degrees  R  using  Equation  1.  Finally, 
the  wall  temperature  on  a  HRSI  tile  was  varied  in  time  from  530  degrees 


to  a  radiation  equilibrium  temperature  where 

•  25 

T'eq  =/  q'\ 

\ae  / 

and  a  is  the  Stefan-Boltzmann  constant  and  e  is  the  tile  emissivity. 

The  above  conditions  may  not  be  an  accurate  calculation  of  the 
small  region  at  the  leading  edge  (12)  and  would  indicate  lower  heating 
near  the  knife  blade  leading  edge  than  the  exact  answer.  The  heating 
near  the  leading  edge  would  also  be  affected  by  a  shock-boundary  layer 
interaction  at  the  lower  deflection  angles.  The  full  or  parabolized 
Navier-Stokes  equations  would  have  to  be  solved  at  each  deflection 
angle  in  order  to  resolve  the  interaction. 

Edge  Boundary  Conditions.  The  boundary  conditions  at  the  outer 
edge  were  approximated  by  oblique  shock  theory.  Oblique  shock  theory 
assumes  inviscid  flow  in  a  perfect  gas.  The  only  pressure  gradients  in 
the  flow  are  assumed  to  be  across  the  oblique  shock.  Therefore,  the 
edge  conditions  were  based  on  the  conditions  behind  the  oblique  shock. 

If  the  plate  is  stationary,  the  edge  conditions  do  not  change.  As  the 
plate  is  pitched,  the  edge  conditions  are  recalculated  for  each  new 
angle  using  Eqs  33  to  37.  Instead  of  solving  the  equations  at  each 
angle,  the  program  interpolates  from  Table  II  which  is  the  result  of  the 
equations  at  each  deflection  angle. 

Initial  Conditions.  To  begin  the  calculations,  an  initial  U  and  V 
velocity  and  total  enthalpy  must  be  known.  The  V  velocity  was  assumed 
to  be  zero  everywhere,  inside  the  boundary  layer,  the  U  velocity  is 
calculated  using  the  cubic  given  as 
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(72) 


5  =  5.2  ST™ 

r^Re1 


for  incompressible  flow  and  y  is  in  the  physical  plane  and  not  the 
transformed  plane.  This  equation  is  used  because  it  gives  a  zero  second 
derivative  at  the  wall  and  gives  a  reasonably  good  velocity  profile 
through  the  boundary  layer  (10).  The  same  cubic  equation  is  used  to 
determine  the  total  enthalpy  profile,  but  it  is  based  on  the 
thermal  boundary  layer  thickness  6T.  The  equation  for  enthalpy  is 


where  ST  =  s/Pr. 

Outside  the  boundary  layer,  U  and  H  are  equal  to  the  edge 


conditions. 


Table  II 


Oblique  Shock  Table 


e  degrees 

Mach 

Velocity 

Temp 

Pressure 

3 

12.203376 

4843.7962 

65.62 

1.3063 

4 

11.541259 

4834.5371 

73.08 

1.7262 

5 

10.886080 

4823.7297 

81.78 

2.2334 

6 

10.245773 

4811.1933 

91.84 

2.8322 

7 

9.630664 

4796.8408 

103.33 

3.5250 

8 

9.046685 

4780.5741 

116.30 

4.3135 

9 

8.498054 

4762.3358 

130.80 

5.1985 

L0 

7.986783 

4742.0933 

146.83 

6.1801 

LI 

7.512550 

4719.8062 

164.39 

7.2581 

L2 

7.0741 

4695.4522 

183.49 

8.4319 

L3 

6.6693 

4669.00 

204.12 

9.7010 

L4 

6.2690 

4639.60 

227 .99 

11.060 

IV.  Results 


To  investigate  the  heat  transfer  on  a  wedge,  different  cases  were 
run  using  the  unsteady  boundary  layer  code  developed  in  this  project. 
The  program  is  verified  using  incompressible  and  compressible  cases. 

Both  cases  are  steady  state  solutions,  each  having  an  isothermal  wall. 
The  incompressible  solution  is  compared  to  the  Blasius  solution  (13:142) 
and  the  compressible  solution  is  compared  to  approximate  isothermal 
Eckert  theory  and  another  numerical  solution  at  three  angles,  3  degrees, 
10  degrees,  and  14  degrees.  After  the  cases  are  verified,  steady  state 
solutions  with  a  nonisothermal  wall  are  run  and  compared  to  theory  and 
experimental  data.  In  addition  to  the  steady  state  solutions,  transient 
cases  are  investigated.  The  transient  cases  change  the  wall  temperature 
of  the  tile  as  a  function  of  time  at  a  constant  of  14  degrees.  Equation 
(1)  is  used  with  a  time  constant  of  10  seconds  to  change  the  wall 
temperature;  in  another  case,  the  change  depends  on  both  the  heat  rate 
and  the  time  constant.  Finally,  the  plate  is  pitched  with  an  isothermal 
wall  in  one  case  and  a  nonisothermal  wall  in  another  case  from  14 
degrees  to  3  degrees  and  back  to  14  degrees. 

Program  Verification 

In  order  to  use  the  program,  THESIS,  developed  in  this  project,  the 
program  must  be  verified.  The  first  step  towards  this  goal  is  running 
an  incompressible  case.  In  this  case,  Mg  is  equal  to  .05  which  is  also 
equal  to  with  a  wall  temperature  of  530  degrees  R.  Since  the 
density  is  essentially  constant,  the  Dorodnitsyn  transformation  is 
negated  and  only  the  incompressible  boundary  layer  equations  are 
verified.  In  Figure  6,  the  numerical  solution  is  compared  to  the 
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530°  R 


Figure  6  Comparison  of  Blasius  Solution  with  the 
Numerical  Solution 


3lasius  solution  (13:42).  The  velocity  ratio,  u/u^,  is  .99  at  equal 
to  5.64  for  the  Blasius  solution  and  6.14  for  the  numerical  solution. 

The  two  curves  compare  well  with  the  max  error  being  five  percent. 

To  verify  the  Dorodnitsyn  transformed  equations,  a  compressible 
case  is  solved  using  the  thin  skin  wind  tunnel  test  as  a  model.  This 
case  is  an  isothermal  wall  at  530  degrees  R  and  I'J,  of  14.237.  The  TQ 
and  PQ  are  2020  degrees  R  and  1600  psia  respectively.  These  conditions 
are  used  in  each  case  for  angles  of  14  degrees,  10  degrees  and  3 
degrees,  and  the  boundary  layer  edge  conditions  are  adjusted  according 
to  oblique  shock  theory. 

Due  to  compressibility  effects,  the  grid  used  in  the  incompressible 
case  must  be  adjusted  to  get  enough  grid  points  in  the  boundary  layer  to 

take  into  account  the  temperature  inversion  and  get  an  accurate 

solution.  This  grid  is  adjusted  by  multiplying  the  Y  location  by  a  grid 
factor.  The  smaller  the  grid  factor,  the  more  points  inside  the 
boundary  layer.  Table  III  shows  the  relationship  of  the  grid  factor  to 
the  h/href  for  the  three  angles.  A  grid  factor  of  1.0  is  used  for  this 
work,  realizing  that  this  is  still  a  relatively  coarse  grid. 

Once  the  grid  is  modified,  a  solution  for  the  compressible  case  can 

be  determined  but  the  result  will  have  a  large  error  due  to  the  large 

gradients  at  the  leading  edge  and  the  coarse  grid.  A  unique  way  of 
taking  care  of  this  leading  edge  problem  is  to  assume  a  zero  gradient  at 
the  leading  edge.  This  implies  that  U(l,j)  is  equal  to  U(2,j),  H (1 , j) 
is  equal  to  H(2,j),  and  V(l,j)  is  equal  to  V(2,j).  This  assumption 
could  be  used  in  other  problems  that  do  not  have  similarity  solutions 
because  the  dx  could  be  made  sufficiently  small  to  make  the  assumption 


Table  lit 


Grid  Factor  Effects  on  h/hre£ 


grid  factor 


3  degrees 

.1 

1.4484 

.5 

1.44375 

1.0 

1.42639 

3.0 

1.28286 

5.0 

1.09215 

approximate 

isothermal 

theory 

1.65 

other 

numerical 

solution 

1.55 

h/\ef 


degrees 

14  degrees 

3.1334 

4.1709 

3.1342 

4.1729 

3.1262 

4.16815 

3.0467 

4.1042 

2.9089 

3.99679 

3.55 

4.6552 

3.35  4.45 


valid.  When  this  assumption  is  used,  the  results  for  the  isothermal 
compressible  case  were  constant  as  the  approximate  isothermal  Eckert 
theory  predicted. 

Upon  the  resolution  of  the  grid  and  the  leading  edge  error,  the 
steady  state  solution  for  the  three  angles  was  calculated.  To  avoid  any 
start  up  problems,  the  14  degrees  steady  state  solution  was  used  as  the 
10  degrees  initial  conditions  and  the  10  degrees  steady  state  solution 
was  used  as  the  3  degrees  initial  conditions.  The  results  are  compared 
in  Figure  7  to  the  approximate  analytical  solution,  to  another  steady 
numerical  boundary  layer  solution  (2),  and  to  in  skin  wedge  data  from 
the  wind  tunnel.  As  is  shown,  the  h/href  results  from  THESIS  are  lower 
than  both  the  approximate  theory  by  approximately  10  percent  and  the 
numerical  solution  by  about  6  percent  at  14  degrees  and  10  degrees  and 
lower  at  3  degrees.  An  explanation  for  this  difference  goes  back  to 
Table  III,  which  shows  that  the  solution  depends  on  the  grid  and 
ultimately  the  grid  factor.  Six  times  as  many  grid  points  were  used  in 
the  other  steady  boundary  layer  solution.  The  thin  skin  data  varies 
more  than  the  difference  in  the  predicted  solutions,  which  can  be  for 
several  reasons.  Variation  in  the  thin  skin  heating  data  along  the 
panel  can  be  due  to  small  nonisothermal  wall  effects  and  thermal 
expansion  of  the  thin  skin  panel.  In  addition,  a  jet  of  air  from  the 
tunnel  start  up  before  model  injection  was  increasing  the  temperature  on 
the  downstream  part  of  all  panels,  and  caused  lower  heating  at  the  rear 
of  the  panel  (7:5).  Some  variation  may  also  be  due  to  shock  interaction 
with  the  viscous  boundary  layer. 
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Figure  7  Thin  Skin  Wedge  Heating  Comparison  with 
Predictions  for  isothermal  Wall 


Even  though  the  numerical  solutions  from  TOESIS  are  grid  dependent, 
the  results  compare  well  with  both  the  approximate  isothermal  Eckert 
theory  and  another  steady  numerical  solution.  This  implies  that  the 
Dorodnitsyn  transformed  boundary  layer  equations  are  correct. 
Nonisothermal  Cases 

To  begin  the  investigation  of  nonisothermal  wall  effects,  some 
steady  state  solutions  are  calculated  and  compared  to  approximate 
analytic  solutions  and  HRSI  and  FRSI  wind  tunnel  data.  The  free  stream 
conditions  are  the  same  as  the  isothermal  case  with  a  ^  of  14.237,  PQ 
of  1600  psia,  and  TQ  of  2020  degrees  R,  but  in  these  cases,  a 
temperature  step  is  given  at  a  streamwise  location  of  .46667.  In  order 
to  compare  wind  tunnel  data  to  the  numerical  results  for  the  HRSI  panel, 
a  temperature  step  of  270  degrees  R  is  given,  making  the  wall  temperature 
jump  from  530  degrees  R  to  800  degrees  R.  The  solution  of  the  boundary 
layer  equations  at  14  degrees  deflection  is  shown  in  Figure  8.  (Versus 
streamwise  distance.)  The  approximate  nonisothermal  solution  given  by 
equation  32  demonstrates  the  dramatic  decrease  in  the  heat  transfer 
coefficient  just  downstream  of  the  temperature  jump,  and  the  subsequent 
slow  recovery.  The  numerical  boundary  layer  solution  shows  the  same 
trend,  but  at  a  slightly  lower  heating  magnitude  similar  to  other 
predictions  (7:6).  The  only  available  surface  thermocouple  data  show  an 
even  lower  magnitude  and  slower  recovery  of  the  boundary  layer  (7:6). 

For  completeness,  solutions  of  the  boundary  layer  equations  at  10  and  3 
degree  deflection  angles  and  the  same  temperature  step  are  shown  in 


Figures  9  and  10  respectively.  No  wind  tunnel  results  are  available  for 
these  two  solutions.  At  3  and  10  degrees,  the  heating  magnitude  is 
slightly  lower  than  the  approximate  nonisothermal  solution. 

To  compare  with  FRSI  wind  tunnel  data,  a  steady  state  nonisothermal 
wall  case  at  a  14  degree  deflection  was  run  with  a  temperature  step  of 
170  degrees  R  making  the  wall  temperature  jump  from  530  degrees  R  to  700 
degrees  R.  The  results,  shown  in  Figure  11,  show  the  same  dramatic 
decrease  in  heat  transfer  coefficient  downstream  of  the  jump  for  both 
the  approximate  nonisothermal  solution  and  the  THESIS  numerical 
solution,  but  are  still  higher  than  wind  tunnel  results  except  at  the 
end  of  the  plate  which  is  lower.  The  boundary  layer  appears  not  to 
recover  as  fast  as  the  predictions  near  the  interface,  but  does  recover 
downstream  reasonably  close  to  predicted  heating  (7:6). 

In  all  of  the  steady  state  nonisothermal  cases,  the  numerical 
results  compared  well  with  the  approximate  nonisothermal  solution  but 
were  consistently  higher  than  the  wind  tunnel  data. 

Transient  Cases 

Understanding  how  the  transient  effects  influence  the  heating  rate 
is  important  for  interpreting  the  overall  results.  To  investigate  these 
effects,  the  wall  temperature  is  increased  as  a  function  of  time  by 
equation  1  from  530  degrees  R  to  800  degrees  R.  The  deflection  angle  is 
held  constant  at  14  degrees.  The  14  degree  steady  state  solution  for 
the  isothermal  wall  is  used  as  the  initial  condition  and  the  freestream 
conditions  are  the  same  as  the  steady  state  nonisothermal  wall  case. 

The  physical  time  step  is  .1  second  and  the  convergence  criteria  are 
specified  as  1.0  e-6.  In  Equation  1,  the  temperature  is  dependent  on 
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Figure  10 


t^=  14.237 

Tw  -  530°  R 

T,,  *  800°  R 

w2 

•  •  •  • 

Approximate  Analytical 

Theory 

THESIS  Numerical  Solution 

i 

— - 

'<  / 

!  / 

T 

0.4  0.5  0.6  0.7  0.8  0 

X'/L 


Non isothermal  Wall  Solution  at  3  Degrees  with 
a  HRSI  Wall  Temperature  of  800  Degrees  R 


both  time  and  a  time  constant.  A  time  constant  of  30  seconds  is  used  to 
simulate  the  conduction  into  the  FRSI  material.  As  the  solution  marches 
in  time  the  wall  temperature  is  changed.  The  solution,  shown  in  Figure 
12,  converges  as  time  increases. 

To  determine  if  a  lag  in  the  h/href  occurred  due  to  the  change  in 
wall  temperatures,  some  additional  steady  state  isothermal  wall  problems 
were  run.  At  a  streamwise  location  of  .5,  the  wall  temperature  at  10 
seconds  is  recorded  to  be  635.5.  This  wall  temperature  was  input  for 
an  isothermal  case  to  determine  the  isothermal  heat  transfer  coefficient 
ratio.  The  isothermal  h/href  ratio  is  3.26207  and  the  transient  h/href 
ratio  is  3.2578  which  indicates  that  a  lag  exists  with  this  large  At  of 
.05. 

To  make  this  transient  case  more  realistic,  the  wall  temperature  is 
coupled  to  the  heat  rate  into  the  wall  and,  as  time  increases,  an 
equilibrium  wall  temperature  will  be  reached.  The  equation  used  to 
calculate  the  wall  temperature  is 

T'w  =  T'w  e-At'/RCI  +  (l-e-At'/RC,)T,eq  (74) 

where 

T'eq  =^q'j1/4  (75) 

where  e  is  .85,  a  is  the  Stephen  Boltzman  constant  4.761  X  10“  ,  and  dt 

is  .05  seconds.  A  time  constant  of  10  seconds  is  used  in  this  case  to 
simulate  the  lag  due  to  conduction  into  a  HRSI  panel.  An  assumption  is 
made  here  that  the  heat  rate  into  the  panel  equals  the  radiation  heat 
rate.  This  is  a  good  assumption  for  shuttle  flight  but  not  for  the  wind 
tunnel;  therefore,  the  equilibrium  wall  temperature  in  this  case  can  be 
higher  than  in  the  wind  tunnel.  The  solution  of  h/href  at  5,  10,  15, 


and  20  seconds  and  steady  state  versus  streamwise  distance  are  shovm  in 
Figure  13.  As  referenced  earlier,  the  equilibrium  temperature  is 
approximately  900  degrees  R  at  the  end  of  the  plate,  which  is  higher 
than  the  800  degrees  R  in  the  wind  tunnel.  This  solution  is  used  for 
the  initial  conditions  of  the  nonisothermal  pitching  plate  case. 

Pitching  Plate  Cases 

Pitching  the  plate  with  the  wall  temperature  as  a  function  of  both 
the  heat  rate  and  time  is  the  most  realistic  case.  A  pitching  plate 
with  an  isothermal  wall  is  examined.  Using  an  isothermal  pitching 
plate,  the  error  due  to  the  first  order  time  difference  is  investigated. 
The  wall  temperature  is  constant  at  530  degrees  R  and  the  freestreain 
conditions  are  the  same  as  the  isothermal  cases.  The  plate  is  pitched 
at  2  degrees  per  second  and  the  edge  conditions  are  changed  accordingly 
with  interpolation  of  oblique  shock  tables.  The  first  At  is  .1  seconds 
with  a  maximum  of  5  iterations  to  converge,  then  the  At  is  halved  to  .05 

with  a  maximum  of  5  iterations  to  converge.  Finally,  a  a  t  of  .025  is 

run  with  maximum  iteration  limits  of  5  and  10. 

The  h/hre£  ratios  are  compared  at  Time  steps  of  .1,  .05,  and  .025 
at  an  a  of  12  degrees  and  an  X/L  location  of  .5.  The  h/href  is  3.45645 
at  At  of  .1,  3.5351  at  A  t  of  .05  and  3.5887  at  A  t  of  .025  which  is  less 
than  2  percent  error  between  A  t  of  .05  and  .025.  Furthermore,  the  At  of 

.025  had  less  than  1  percent  error  between  5  and  10  iterations.  This 

implies  that  the  unsteady  terms  could  have  been  ignored  without  picking 
up  more  than  1%  error  if  a  small  enough  time  step  is  used  and  the 
solution  is  a  steady  state  solution  at  each  angle.  With  this  in  mind,  a 
r  function  (8:150)  could  have  been  used  to  solve  for  the  heat  rate  by 
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q  =  .332  I^Pr1/3  Re1/2  (  I  nBr4  x*n  ^  +  A) 
x*  x 

where 

8  =  r(4/3n)  r (2/3)  (77) 

r(4/3n  +  2/3) 

Finally,  note  the  h/href  is  lower  at  the  leading  edge  and  recovers  to  a 
constant  at  the  end  of  the  plate.  The  response  is  faster  at  the  leading 
edge  because  the  boundary  layer  is  thinner  and  it  takes  less  time  for 
energy  to  diffuse  across  the  boundary  layer.  It  also  takes  less  time 
near  the  leading  edge  for  energy  to  convect  downstream  from  any  changes 
upstream. 

Pitching  the  plate  with  a  nonisothermal  wall  emulates  pitching  the 
plate  in  the  wind  tunnel  and  the  push  over  pull  up  maneuver  used  by  the 
Space  Shuttle  (7).  The  wall  temperature  is  a  function  of  time  and  heat 
rate  where  Equations  74  to  75  apply.  The  time  constant  is  changed  to  3 
seconds  to  emulate  the  material  response  of  the  HRSI  tile  during  a  small 
maneuver  and  the  converged  solution  of  the  transient  wall  temperature  as 
a  function  of  heat  rate  is  used  as  initial  conditions.  As  before,  the 
plate  is  pitched  at  2  degrees  per  second  with  a  deflection  angle  time 
history  shown  in  Figure  14.  Starting  at  equilibrium,  the  plate  is 
pitched  from  14  degrees  to  3  degrees  and  held  for  5  seconds.  Then  it  is 
pitched  back  up  to  14  degrees  and  held  for  another  5  seconds.  This 
solution  has  some  very  interesting  results.  Figures  15  and  16  are  plots 
of  the  h/href  and  wall  temperature,  respectively,  versus  streamwise 
location  for  14,  9,  and  3  degrees  as  the  plate  is  pitched  up  and  3,  9, 
and  14  degrees  as  it  is  pitched  down.  Both  Figures  show  a  lag  in  the 
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wall  response,  particularly  at  the  end  of  the  plate,  which  is  reasonable. 
A  small  transient  effect  can  be  seen  in  Figure  15  near  the  leading  edge. 
In  addition.  Figure  15  shows  the  same  trend  as  previously  discussed  in 
the  nonisothermal  wall  cases.  However,  at  the  interface,  the  h/hre£  are 
about  the  same  for  all  angles  before  they  changed.  Furthermore,  in  many 
of  the  previous  cases,  a  temperature  step  was  assumed  and  Figure  16  is  a 
good  illustration  that  the  assumption  was  an  acceptable  one. 

In  addition  to  these  two  plots,  the  h/href  and  wall  temperature 
time  histories  and  the  h/hre£  versus  the  deflection  angle  are  plotted 
for  five  streamwise  locations,  .016667,  .466667,  .483333,  .566667,  and 
.983333.  The  first  point,  .016667,  is  right  at  the  leading  edge  of  the 
plate  and  Figures  17,  18,  and  19  are  the  above  mentioned  plots  for  this 
location.  Figure  17  shows  that  the  h/hre£  decreases  as  the  plate  is 
pitched  to  3  degrees  and  increases  as  it  pitches  back  down  to  14  degrees 
due  to  the  change  in  the  heat  rate.  At  the  5  second  hold  at  both  3 
degrees  and  14  degrees  the  solution  is  at  steady  state  almost 
immediately,  but  the  unsteady  terms  in  the  equations  do  affect  the 
results  at  about  10.5  and  16  seconds.  Figure  18  shows  the  temperature 
is  a  constant  530  degrees  R  to  simulate  the  steel  wind  tunnel  model 
throughout  the  plate  pitching,  which  implies  that  there  are  no  wall 
effects  at  this  location.  The  Vhre£  versus  the  deflection  angle. 

Figure  19,  indicates  a  small  hysteresis.  This  could  have  numerical 
error  due  to  the  large  gradients  in  the  boundary  layer  at  the  leading 
edge  and  the  coarse  grid. 


x’/L 


Figure  15  Comparison  of  Nonisothermal  Pitching  Plate  Solutions 
for  14  Degrees,  9  Degrees,  and  3  Degrees 
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At  the  interface  of  the  panel,  .46667  location.  Figures  20,  21,  and 
22  look  very  similar  to  the  .016667  location.  Since  conduction  through 
the  plate  is  not  accounted  for  in  this  project,  the  wall  temperature 
remains  at  530  degrees  R,  Figure  21,  even  though  it  is  next  to  a  higher 
or  lower  temperature  depending  on  the  direction  of  the  pitch.  Further¬ 
more,  the  hysteresis.  Figure  22,  is  less  than  at  the  leading  edge  but 
that  is  to  be  expected  since  the  boundary  layer  gradients  have 
decreased.  Therefore  this  hysteresis  is  probably  caused  by  the  error 
due  to  the  coarse  grid  and  some  unsteady  effects. 

For  the  first  location  after  the  interface  of  the  panel,  .483333, 
the  h/hre£  changes  very  little  as  the  plate  is  pitched.  Figure  23.  The 
equilibrium  h/hre£  is  much  lower  than  the  previous  location  because  the 
wall  temperature  is  about  200  degrees  higher,  Figure  24.  Hysteresis 
also  exists  at  this  location,  Fig  25. 

Further  down  the  plate,  at  location  .56667,  the  results  are  much 
more  dramatic,  Figures  26,  27,  and  28.  As  the  plate  is  pitched  the 
h/href  decreases.  When  the  deflection  angle  first  reaches  3  degrees, 
the  wall  temperature  is  higher  than  the  equilibrium  wall  temperature 
because  of  a  lag.  While  the  plate  remains  at  3  degrees  for  5  seconds, 
the  wall  temperature  decreases  and,  in  tum,  the  h/hre£  increases  until 
equilibrium  is  reached.  The  five  second  hold  was  not  quite  long  enough 
to  reach  equilibrium.  The  same  phenomenon  happens  as  the  plate  is 
pitched  back  to  14  degrees,  only  in  reverse.  When  the  plate  reaches  14 
degrees  the  wall  temperature  is  lower  and  then  rises  to  equilibrium 
during  the  following  five  seconds.  As  the  temperature  rises,  the  h/hrejr 
decreases.  Some  unsteady  effects,  noticed  in  the  solution  of  the 
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Figure  22  h/href  versus  Deflection  Angle  for  x/L  Location 
of  .466667 
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leading  edge  o£  the  plate,  may  be  incorporated  in  this  solution  but  it 
is  impossible  to  separate  them.  Finally,  the  results  show  a  large 
hysteresis  at  this  location;  which  indicates  a  large  time  lag.  This  lag 
is  mainly  due  to  the  large  temperature  gradients  at  this  location  and 
possibly  some  error  due  to  the  coarse  grid. 

Finally,  near  the  end  of  the  plate  at  a  streamwise  location  of 
.983333,  Figures  29,  30,  and  31,  the  solution  has  relatively  recovered. 
Figure  29  looks  much  like  Figure  17  only  at  a  lower  heating  rate.  Some 
nonisothermal  wall  effects  on  Figure  29  are  still  noticeable  at  the  3 
degrees  and  14  degrees  five  second  hold  and  are  caused  for  the  same 
reasons  previously  described.  When  h/hre£  is  plotted  versus  the 
deflection  angle,  Figure  31,  a  lag  due  to  nonisothermal  effects  still 
exists  but  less  than  at  the  previous  location. 

In  conclusion,  this  pitching  plate  case  with  the  wall  temperature, 
as  a  function  of  the  heat  rate,  does  show  a  time  lag  due  to  nonisotherm- 
al  wall  effects.  As  the  plate  pitches  the  wall  temperature  lags  and 
when  the  plate  is  held  at  a  constant  angle,  the  wall  temperature  changes 
until  it  reaches  equilibrium  which  also  causes  the  h/hreE  to  change. 
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Figure  31 

h/href  versus  Deflection  Angle 
of  .983333 

for  x/L  Location 

V.  Conclusions  and  Recommendations 

Conclusions 

The  significant  findings  of  this  research  can  be  put  into  two 
categories.  One,  the  unique  features  of  the  code  written  to  investigate 
nonisothermal  wall  effects,  and  two,  the  actual  heat  transfer  results. 

The  code  developed  in  this  project  solves  the  unsteady  compressible 
boundary  layer  equations.  In  code  development,  the  boundary  layer 
equations  are  transformed  using  the  Dorodnitsyn  transformation.  This  is 
a  unique  feature  because  the  equations  can  be  solved  much  like  the 
incompressible  equations.  These  equations  are  then  transformed  to  the 
computational  plane  and  finite  differenced.  A  three  point  upwind 
differencing  scheme  is  used  for  the  convective  derivatives,  a  second 
order  central  difference  is  used  for  viscous  terms,  and  the  time  terms 
use  a  first  order  backward  difference.  Finally,  the  momentum  and  energy 
equations  are  solved  using  optimized  successive  over-relaxation,  and  the 
continuity  equation  is  solved  for  the  V  velocity  by  integration. 

Another  unique  feature  of  this  code  is  a  generalized  grid  (as  long 
as  Xn  =0).  This  enables  the  program  to  handle  various  configurations. 
Although,  if  another  grid  is  used,  another  method  of  handling  the  Y^ 
will  have  to  be  used  because,  at  the  present  time,  the  analytical 
solution  of  Yj-  based  on  the  Blasius  grid  is  used. 

The  grid  used  in  the  project  is  based  on  an  optimized  Blasius  grid. 
When  a  compressible  case  is  run,  the  grid  must  be  adjusted  by  a  grid 
factor  to  obtain  enough  points  in  the  thermal  boundary  layer.  Even  with 
the  adjustment,  there  are  only  61  points  in  the  streamwise  direction  and 


30  points  in  the  Y  direction  that  go  out  a  distance  of  10  s.  This  code 
has  significantly  fewer  grid  points  than  other  schemes  which  will  enable 
the  code  to  run  faster  and  less  costly. 

The  boundary  conditions  also  played  an  important  part  in  this 
research.  For  the  most  complicated  case,  pitching  a  plate  with  a 
transient  nonisothermal  wall,  both  the  edge  conditions  of  the  boundary 
layer  and  the  wall  temperature  would  change  as  the  plate  is  pitched. 

A  final  feature  of  this  code  is  the  way  the  leading  edge  start  up 
is  handled.  When  the  code  was  first  used  for  an  isothermal  case,  error 
due  to  the  leading  edge  was  seen.  A  zero  gradient  is  assumed  at  the 
leading  edge  and  this  greatly  reduced  the  error.  This  similarity  at  the 
leading  edge  should  be  applicable  to  other  problems  if  a  small  enough 
step  is  taken  in  the  streamwise  direction. 

The  code  was  used  to  solve  for  many  cases,  such  as,  an  incompress¬ 
ible  case,  steady  state  isothermal  and  nonisothermal  cases,  transient 
wall  temperature  cases,  and  pitching  plate  cases,  and  the  results  are 
significant.  The  incompressible  case  and  the  isothermal  cases  were 
used  as  a  program  verification.  The  incompressible  case  compared  well 
with  the  Blasius  solution  and  the  isothermal  cases  well  with  the 
approximate  isothermal  Eckert  theory.  From  the  isothermal  cases  it  can 
be  seen  that  the  code  solutions  are  very  grid  dependent.  The 
nonisothermal  cases  followed  the  same  trend  as  the  approximate 
nonisothermal  solution  given  by  Eqn  32.  It  demonstrates  the  same 
dramatic  decrease  in  the  heat  transfer  coefficient  just  downstream  of 
the  temperature  jump  and  the  subsequent  slow  recovery.  The  transient 
nonisothermal  wall  displayed  the  same  trend.  In  addition,  the  transient 


cases  both  converged  as  time  is  marched.  No  lag  due  to  change  in  wall 
temperature  was  evident  but  the  time  step  may  have  been  too  big  to  see 
the  small  effect.  An  interesting  outcome  of  the  transient  case  with  the 
wall  temperature  being  a  function  of  the  heat  rate  is  that  the  wall 
temperature  is  higher  than  wind  tunnel  wall  temperature.  Finally,  the 
pitching  plate  cases  were  examined.  The  pitched  plate  with  an 
isothermal  wall  temperature  was  used  to  investigate  the  error  due  to  the 
first  order  time  difference.  Between  .05  and  .025  seconds,  there  was 
less  than  2  percent  change  for  a  X/L  location  of  .5  and  an  a  of  12 
degrees.  The  pitching  plate  with  the  wall  temperature  as  a  function  of 
the  heat  rate  yielded  some  very  significant  results.  A  significant 
hysteresis  is  not  evident  at  the  beginning  of  the  plate  and,  after  the 
interface  a  hysteresis  due  to  nonisothermal  effects  does  exist,  more 
dramatically  close  to  the  interface  than  at  the  end  of  the  plate. 

This  research  has  been  significant  in  understanding  the 
nonisothermal  wall  effects  for  the  FRSI  and  HRSI  used  on  the  Shuttle. 
Other  locations  on  the  Shuttle  can  be  identified  where  nonisothermal 
wall  effects  will  impact.  Even  though  this  project  was  based  on  a  Space 
Shuttle  application,  the  code  could  and  should  be  used  in  other 
applications. 

Recommendations 

The  following  is  a  list  of  recommendations  set  in  an  order  of 
priority,  even  though,  many  of  the  items  could  be  investigated 
independently. 

1)  The  grid  should  be  optimized  for  a  nonisothermal  wall. 


2)  A  turbulence  model  should  be  added  so  the  full  boundary 
layer  regime  can  be  investigated. 

3)  Modify  the  equations  for  axisymmetric  problems  and, 
along  the  same  lines,  make  the  equations  three  dimen¬ 
sional. 

4)  Add  the  pressure  gradient  term  to  the  momentum  equation 
which  was  assumed  to  be  zero  in  these  applications. 

5)  Using  the  same  Dorodnitsyn  transformation  and  the  simil¬ 
arity  solution  at  the  leading  edge,  investigate  shock¬ 
fitting  and  shock  capturing  schemes  for  the  parabolized 
Navier-Stokes  Equations. 

6)  Investigate  different  finite  difference  solvers  such  as 
the  tridiagonal  solver. 

7)  Couple  this  program  with  a  conduction  model  to  allow 
conduction  through  the  wedge. 
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Dorodnitsyn's  Transformation  of  the  boundary  Layer  Equations 


This  project  is  based  on  the  solution  of  the  boundary  layer 
equations.  The  Dorodnitsyn  transformation  allowed  the  compressible 
unsteady  boundary  layer  equations  to  be  solved  similar  to  the 
incompressible  unsteady  equations.  The  details  of  the  transformation 
are  included  here  for  reference. 

The  unsteady  compressible  boundary  layer  equations  are 
nondimensionalized  by 
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in  nonconservative  form.  The  Reynolds  number  (Re)  and  Prandtl  number 
(Pr)  are  defined  as 
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is  used  to  transform  the  equations,  A-2  to  A-4,  into  a  form  similar  to 


the  incompressible  equations.  Using  the  chain  rule,  jhj  can  be  defined 
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where  _3Y  =  p  and  _3t'  = 
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Substituting  these  quantities  into  the  continuity  equation  gives 
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equation  becoming  like  the  incompressible  form 
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The  same  procedure  is  used  for  the  momentum  equation.  Beginning 
with  the  nonconservative  momentum  equation,  Eqn  A-3,  and  again 
substituting  and  canceling  terms,  the  equation  becomes 
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